Action

Building on the momentum and ideas of our last meeting, the following pipeline was attempted to produce the wildfire fuel mapping outputs described in the NRC grant “High-Resolution Mapping”:

Using the new cffdrs R-package (@wang2017cffdrs; @van1985equations; @van1987development), we developed rasters of forest fuel moisture code and wildfire weather indices (Table 1) that were used to fit the stand-adjusted fine-fuel model (@wotton2007stand) applied to vegetation rasters classified according to 16 fuel classes of the BC forest fuel typing algorithm ( @perrakis2018british). REeult rasters were then used to fit the Canadian Forest Fire Behaviour Prediction model from which we drafted raster maps representing Head Fire Index (HFI) and Fire Intensity maps (FI) (Figure 1) for the Okanagan Watershed Basin for the day of June 30th 2021.

Figure 1: Tentative Workflow

0.1 Import Site Data

We applied the Okanagan Watershed boundary (FWA ID:212) as our area of interest, which was downloaded from the BC Geographic Warehouse. We tested potential base mapping services from Google Cloud API and the ggmap package with a free personal account key token. Four basemaps were derived at a zoom setting of 9 at the latitude and longitude location of (-119, 50.0). Likely that Google Cloud services are working off EPSG 3857, though this needs confirming. With the default ‘statellite’ basemap, we used Lobo’s script (2014) to transform the RGB raster to a matrix object and a SpatRaster before applying different mapping styles. Note that chunk outputs require substantial system memory and may cause crashes.

register_google(key = 'AIzaSyDbjpyUjJp3VInMOqebU9yp2zff5hA-zBM')
gmap1 = ggmap(get_map(location = c(-119.7, 50.0), maptype = "satellite", source = "google", zoom = 9))
gmap2 = ggmap(get_map(location = c(-119.7, 50.0), maptype = "toner-lite", zoom = 9))
gmap3 = ggmap(get_map(location = c(-119.7, 50.0), maptype = "toner-background", zoom = 9))
gmap4 = ggmap(get_map(location = c(-119.7, 50.0), mqptype="terrain-labels", zoom = 9))

0.2 Google Cloud Api & GUI Basemap Pipelines

mgmap <- as.matrix(gmap)
vgmap <- as.vector(mgmap)
vgmaprgb <- col2rgb(vgmap)
gmapr <- matrix(vgmaprgb[1, ], ncol = ncol(mgmap), nrow = nrow(mgmap))
gmapg <- matrix(vgmaprgb[2, ], ncol = ncol(mgmap), nrow = nrow(mgmap))
gmapb <- matrix(vgmaprgb[3, ], ncol = ncol(mgmap), nrow = nrow(mgmap))
rgmaprgb <- brick(raster(gmapr), raster(gmapg), raster(gmapb))
rm(gmapr, gmapg, gmapb)
rgmaprgbGM <- rgmaprgb
raster::projection(rgmaprgbGM) <- CRS("+init=epsg:3857")
extent(rgmaprgbGM) <- unlist(attr(gmap, which = "bb"))[c(2, 4, 1, 3)]
unlist(attr(gmap, which = "bb"))[c(2, 4, 1, 3)]
rprobextSpDF <- as(extent(unlist(attr(gmap, which = "bb"))[c(2, 4, 1, 3)]), "SpatialPolygons")
raster::projection(rprobextSpDF) <- CRS("+init=epsg:4326")
rprobextGM <- spTransform(rprobextSpDF, CRS("+init=epsg:3857"))
extent(rgmaprgbGM) <- c(rprobextGM@bbox[1, ], rprobextGM@bbox[2, ])
rgmaprgbGM_disksafe = writeRaster(rgmaprgbGM, file = "rgmaprgbGM", format = "GTiff", overwrite = TRUE, datatype = "INT1U")
xcenter <- (extent(rgmaprgbGM)@xmax + extent(rgmaprgbGM)@xmin)/2
ycenter <- (extent(rgmaprgbGM)@ymax + extent(rgmaprgbGM)@ymin)/2
height <- extent(rgmaprgbGM)@xmax - extent(rgmaprgbGM)@xmin
width <- extent(rgmaprgbGM)@ymax - extent(rgmaprgbGM)@ymin

The Okanagan Watershed boundary was extracted from British Columbia Freshwater Atlas Dataset and processed as a simple feature reprojected to the EPSG:3857 coordinate reference system match the Google Cloud Mapping API.

watershed_okanagan = st_read("./Data/watershed-okanagan-QX.shp")
## Reading layer `watershed-okanagan-QX' from data source 
##   `/Volumes/128GB_WORKD/EFI-FIRE/01_Wildfire-Mapping-CFFDRS-2.0_prototype/Data/watershed-okanagan-QX.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 1413784 ymin: 464311.2 xmax: 1515562 ymax: 647056.5
## Projected CRS: NAD83(CSRS) / BC Albers
ggplot(watershed_okanagan) + geom_sf(alpha=0.1) + coord_sf()

0.3 Import Topographical Data

LiDAR data of 3-arc second resolution was acquired using the ‘elevatr’ package and the continuing Mapzen license (@van2001shuttle). DEM data was transformed into a spatRaster and disaggreated from 98m resolution fown to 32m~ resolution. Slope and aspect rasters were calculated using the terra::terrain function with a bilinear interpolation and a rook neighbourhood sampling of 8 adjacent cells. This produced problematic results and the raster processing approach was oapplied using the deprecated slopeAspect functions.

ELV = get_elev_raster(watershed_okanagan, z=9)
#raster::projection(ELV) <- CRS("+init=epsg:3857")
GS = slopeAspect(ELV, filename = "./Data/GS.tif", out='aspect', unit='degrees', neighbors=8, overwrite=TRUE)  
GS = raster("./Data/GS.tif")
GS = rast(GS)
Aspect = slopeAspect(ELV, filename = "./Data/Aspect.tif", out='slope', unit='degrees', neighbors=8, overwrite=TRUE)  
Aspect = raster("./Data/Aspect.tif")
Aspect = rast(Aspect)
ELV = rast(ELV)
ELV = disagg(ELV, fact=3.3)
GS = resample(GS, ELV, method="bilinear")
Aspect = resample(Aspect, ELV, method="bilinear")
Aspect = mask(Aspect, vect(watershed_okanagan))
GS = mask(GS, vect(watershed_okanagan))
ELV = mask(ELV, vect(watershed_okanagan))
plot(GS, main='slope')
plot(ELV, main='elevation')
plot(Aspect, main='aspect')

0.4 Import Climate Data

Climate variables were downloaded as NetCDF files from NASA Power platform and read directly into R as rasters of 1) mean daily temperature at 2m, 2) mean daily precipitation, 3) mean relative humidity, 4) and mean wind speed at 10m. The NASA Power platform supports some very user-friendly API links for static data sources that might be useful for this proposed grant project. REminder tho, some API’s prefer dealing with dataframe inputs so might be worht preparing df pipe while still fresh in the head.

temp = terra::rast("./Data/temp.nc")
prec = terra::rast("./Data/prec.nc")
rh = terra::rast("./Data/rh.nc")
ws = terra::rast("./Data/ws.nc")
temp = mask(temp, vect(watershed_okanagan))
prec = mask(prec, vect(watershed_okanagan))
rh = mask(rh, vect(watershed_okanagan))
ws = mask(ws, vect(watershed_okanagan))
temp = mean(temp)
prec = mean(prec)
rh = mean(rh)
ws = mean(ws)
names(temp) = 'temp'
names(prec) = 'prec'
names(rh) = 'rh'
names(ws) = 'ws'
plot(temp, main='temperature (2m)')
plot(prec, main='precipitation (mm/day)')
plot(rh, main='relative humidity')
plot(ws, main='wind speed (10m)')
#temp = terra::resample(temp, ELV, method="bilinear")
#prec = terra::resample(prec, ELV, method="bilinear")
#rh = terra::resample(rh, ELV, method="bilinear")
#ws = terra::resample(ws, ELV, method="bilinear")

0.5 Derive CFFDRS Wildfire Weather Rasters

Interpolated climate predictors were assembled as a raster stack and inputted into the fwiRaster function. The ’out=“all’” option was selectedin the fwiRaster function which gave us the additional raster outputs for Initial Spread Index (isi), and Build-up Index (bui).

temp = raster::raster(temp)
prec = raster::raster(prec)
rh = raster::raster(rh)
ws = raster::raster(ws)
stack = stack(temp, rh, ws, prec)
names(stack)
## [1] "temp" "rh"   "ws"   "prec"
fwi_outputs = fwiRaster(stack, out = "all")
plot(fwi_outputs)

ffmc = raster(fwi_outputs, layer=5)
dmc = raster(fwi_outputs, layer=6)
dc = raster(fwi_outputs, layer=7)
isi = raster(fwi_outputs, layer=8)
bui = raster(fwi_outputs, layer=9)
fwi = raster(fwi_outputs, layer=10)
dsr = raster(fwi_outputs, layer=11)

1 Import VRI Data

The VRI dataset was downloaded from imapBC as a shapefile.shp and transformed into simple feature for processing oeprations using sf and dplyr functions. For CFFDRS data requirements, we consulted the two package-provied sample datasets ‘test_fwi’ and ‘test_fpb’ presented below:

library(cffdrs)
print(as_tibble(test_fwi), n = 10)
## # A tibble: 48 × 9
##     long   lat    yr   mon   day  temp    rh    ws  prec
##    <int> <int> <int> <int> <int> <dbl> <int> <int> <dbl>
##  1  -100    40  1985     4    13  17      42    25   0  
##  2  -100    40  1985     4    14  20      21    25   2.4
##  3  -100    40  1985     4    15   8.5    40    17   0  
##  4  -100    40  1985     4    16   6.5    25     6   0  
##  5  -100    40  1985     4    17  13      34    24   0  
##  6  -100    40  1985     4    18   6      40    22   0.4
##  7  -100    40  1985     4    19   5.5    52     6   0  
##  8  -100    40  1985     4    20   8.5    46    16   0  
##  9  -100    40  1985     4    21   9.5    54    20   0  
## 10  -100    40  1985     4    22   7      93    14   9  
## # … with 38 more rows
print(as_tibble(test_fbp), n = 10)
## # A tibble: 20 × 24
##       id FuelType   LAT  LONG   ELV  FFMC   BUI    WS    WD    GS    Dj    D0
##    <int> <fct>    <int> <int> <int> <dbl> <int> <dbl> <int> <int> <int> <int>
##  1     1 C-1         55   110    NA  90     130  20       0    15   182    NA
##  2     2 C2          50    90    NA  97     119  20.4     0    75   121    NA
##  3     3 C-3         55   110    NA  95      30  50       0     0   182    NA
##  4     4 C-4         55   105   200  85      82   0      NA    75   182    NA
##  5     5 c5          55   105    NA  88      56   3.4     0    23   152   145
##  6     6 C-6         55   105    NA  94      56  25       0    10   152   132
##  7     7 C-7         50   125    NA  88.8    15  22.1   270    15   152    NA
##  8     8 D-1         45   100    NA  98     100  50     270    35   152    NA
##  9     9 M-1         47    85    NA  90      40  15.5   180    25   182    NA
## 10    10 M-2         63   120   100  97     150  41     180    50   213    NA
## # … with 10 more rows, and 12 more variables: hr <dbl>, PC <int>, PDF <int>,
## #   GFL <dbl>, cc <int>, theta <int>, Accel <int>, Aspect <int>, BUIEff <int>,
## #   CBH <lgl>, CFL <lgl>, ISI <int>

Wotton and Beverly’s model of stand-adjusted fine fuel moisture content requires five predictor variables, two of which were extracted directly from spatial layers of the the VRI dataset including stand-type (‘SPEC_CD_1 ==?’ & ’SPEC_PCT_1 > 0.80)BCLCS_LV_1,

vri2020_sf = st_read("./Data/BCGW_7113060B_1645786298548_3276/VEG_COMP_LYR_R1_POLY/VEG_R1_PLY_polygon.shp")
st_crs(watershed_okanagan) = 3005
vri2020_sf = st_intersection(st_make_valid(vri2020_sf), watershed_okanagan)
str(vri2020_sf)

fuel_attributes = vri2020_sf %>%
  dplyr::select(BCLCS_LV_1, BCLCS_LV_2, BCLCS_LV_3, BCLCS_LV_4, BCLCS_LV_5, SHRB_HT, SHRB_CC, HERB_TYPE, HERB_COVER, HERB_PCT, NON_VEG_1, BEC_ZONE, BEC_SZONE, N_LOG_DATE, DEAD_PCT, HRVSTDT, SITE_INDEX, SPEC_CD_1, SPEC_CD_2, SPEC_PCT_1, SPEC_PCT_2, PROJ_HT_1, SPEC_AGE_1, DEAD_PCT, STEM_HA_CD, DEAD_STEMS, NVEG_COV_1)

library(dplyr)
Wotton_fuel_N = 0
Wotton_fuel_decid = 1
Wotton_fuel_Df = 2
Wotton_fuel_MW = 3
Wotton_fuel_PI = 4
Wotton_fuel_SP = 5

Wotton_fuel_class = vri2020_sf%>% 
  mutate(fuel_type = case_when((BCLCS_LV_1 == "N") ~ 0,
  (BCLCS_LV_1 == "V" & BCLCS_LV_4 = "TB") ~ 1,
  (BCLCS_LV_1 == "V" & SPEC_CD_1 == "FD" | SPEC_CD_1 == "FDC" | SPEC_CD_1 == "FDI") ~ 2,
  (BCLCS_LV_1 == "V" & SPEC_PCT_1 <= 80) ~ 3,
  (BCLCS_LV_1 == "V" & SPEC_CD_1 == "PA" | SPEC_CD_1 == "PL" | SPEC_CD_1 == "PLC" | SPEC_CD_1 == "PY") ~ 4,
  (BCLCS_LV_1 == "V" & BCLCS_LV_5 =="SP" | SPEC_CD_1 == "SB" | SPEC_CD_1 == "SX" | SPEC_CD_1 == "SW" | SPEC_CD_1 == "S" | BEC_ZONE == "BWBS" | BEC_ZONE == "SWB") ~ 5,
  TRUE ~ 0))


density = vri_sf["LIVE_STEMS"] %>% mutate(LIVE_STEMS = as.numeric(LIVE_STEMS))
density = rename(density, density = LIVE_STEMS)
ggplot(stand) + geom_sf(aes(fill=stand), size = 0.05)
ggplot(density) + geom_sf(aes(fill=density), size = 0.0005) + scale_fill_viridis_c()
vri2020_sf$SPEC_PCT_1
#M1_fuel_mixedwood80 = vri2020_sf$
stand = vri_sf["SPEC_CD_1"] %>% mutate(SPEC_CD_1 = as.factor(SPEC_CD_1))


raster::crs(watershed_bcimap) = "EPSG:3005"
plot(watershed_bcimap)
watershed_bcimap$Export


plot(st_geometry(vri_sf))
plot(st_geometry(cutblk_sf))
fields::stats(vri_sf) 
fields::stats(cutblk_sf) 

fields::stats(vri_sf$HRVSTDT) 
fields::stats(vri_sf$C_I_CODE) 
fields::stats(vri_sf$BEC_ZONE) 
fields::stats(vri_sf$BCLCS_LV_1)

C1_fuel = vri_sf %>% 
  dplyr::filter(BCLCS_LV_2 == 'V', BCLCS_LV_1 == 'V', 
                HRVSTDT > 19970000 | HRVSTDT > 20140000 | HRVSTDT == NA, 
                BEC_ZONE == "BWBS" \| BEC_ZONE == "SWB", 
                SPEC_CD_1 =='S' | SPEC_CD_1 == 'SB' | SPEC_CD_1 == 'SE' | SPEC_CD_1 =='SX')

M1_fuel = vri_sf %>% d

C2_fuel = vri_sf %>% dplyr::filter(BCLCS_LV_2 == 'V', BCLCS_LV_1 == 'V', HRVSTDT \> 19970000, HRVSTDT \> 20140000, BEC_ZONE == "BWBS" \| BEC_ZONE == "SWB", SPEC_CD_1 =='S' \| SPEC_CD_1 == 'SB' \| SPEC_CD_1 == 'SE' \| SPEC_CD_1 =='SX')

C3_fuel = vri_sf %>% dplyr::filter(BCLCS_LV_2 == 'V', BCLCS_LV_1 == 'V', HRVSTDT \> 19970000, HRVSTDT \> 20140000, SPEC_CD_1 =='Pl' \| SPEC_CD_1 == 'Pli' \| SPEC_CD_1 == 'Plc' \| SPEC_CD_1 =='Pj'\| SPEC_CD_1 == 'P' \| SPEC_CD_1 =='SE',)

C4_fuel = vri_sf %>% dplyr::filter(BCLCS_LV_2 == 'V', BCLCS_LV_1 == 'V', HRVSTDT \> 19970000, HRVSTDT \> 20140000, SPEC_CD_1 =='Pl' \| SPEC_CD_1 == 'Pli' \| SPEC_CD_1 == 'Plc' \| SPEC_CD_1 =='Pj'\| SPEC_CD_1 == 'P' \| SPEC_CD_1 =='SE',)

C4_fuel = vri_sf %>% dplyr::filter(BCLCS_LV_2 == 'V', BCLCS_LV_1 == 'V', HRVSTDT \> 19970000, HRVSTDT \> 20140000, SPEC_CD_1 =='Pl' \| SPEC_CD_1 == 'Pli' \| SPEC_CD_1 == 'Plc' \| SPEC_CD_1 =='Pj'\| SPEC_CD_1 == 'P', vri_sf\$PROJ_AGE_1 \< 2) \# Competitor Tools for Fuel Typing in BC 2022

#### *fwiRaster and sdmc calculated based on daily climate records*

#### *gfmc and hffmc calculated based on hourly climate records - key to CFFDRSv2.0*

#### *Start date of fire season calculated with fireSeason*

#### *All outputs generated for once-daily calcuylations for the full fireSeason chronologically using 'batch=TRUE' function*